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Abstract. Massive black holes are key ingredients of the assembly and evolution 
of cosmic structures. Pulsar Timing Arrays (PTAs) currently provide the 
only means to observe gravitational radiation from massive black hole binary 
systems with masses > 10'^ M0. The whole cosmic population produces a signal 
consisting of two components: (i) a stochastic background resulting from the 
incoherent superposition of radiation from the all the sources, and (ii) a handful 
of individually resolvable signals that raise above the background level and 
are produced by sources sufficiently close and/or massive. Considering a wide 
range of massive black hole binary assembly scenarios, we investigate both the 
level and shape of the background and the statistics of resolvable sources. We 
predict a characteristic background amplitude in the interval hc{f = 10~* Hz) ^ 
5 X 10"^'' — 5 X 10~^^, within the detection range of the complete Parkes PTA. On 
average, at least one resolvable source produces timing residuals that integrated 
over the typical time of observation lay in the range ~ 5 — 50 ns.We also 
quantify the capability of PTAs of measuring the parameters of individual sources, 
focusing on the astrophysically more likely monochromatic signals produced by 
binaries in circular orbit. We investigate how the results depend on the number 
and distribution of pulsars in the array, by computing the variance-covariance 
matrix of the parameter measurements. For plausible Square Kilometre Array 
(SKA) observations (100 pulsars uniformly distributed in the sky), and assuming 
a coherent signal-to-noise ratio of 10, the sky position of massive black hole 
binaries can be located within a 40deg^ error box, opening promising prospects 
for detecting a putative electromagnetic counterpart to the gravitational wave 
emission. The planned SKA, can plausibly observe these unique systems, although 
the number of detections is likely to be small. These observations would naturally 
complement on the high-mass end of the black hole distribution function future 
surveys carried out by the Laser Interferometer Space Antenna (LISA). 



1. Introduction 



Massive black hole (MBH) binary (MBHB) systems with masses in the range ~ 
10^—10"'^° Mq are amongst the primary candidate sources of gravitational waves (GWs) 
at ~ nHz - mHz frequencies [1] [2l |3l |4l [5] . The frequency band ~ 10"^ Hz — IHz will 
be probed by the Laser Interferometer Space Antenna {LISA, [6]), a space-borne GW 
laser interferometer being developed by ESA and NASA. The observational window 
10~® Hz — 10~® Hz is already accessible with Pulsar Timing Arrays (PTAs; e.g. the 
Parkes radio-telescope, [7]). PTAs exploit the effect of GWs on the propagation of 
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radio signals from a pulsar to the Earth [SI HJ [TU] , producing a characteristic signature 
in the time of arrival (TOA) of radio pulses. The timing residuals of the fit of the 
actual TOAs of the pulses to the TOAs predicted by a given pulse emission and 
propagation model carry the physical information about unmodelled effects, including 
GWs [TTl [12] . The complete Parkes PTA [7] , the European Pulsar Timing Array 
[13], and NanoGrav [M] are expected to improve considerably on the capabilities of 
these surveys, and the planned Square Kilometer Array (SKA; www.skatelescope.org) 
will produce a major leap in sensitivity. 

Popular scenarios of MBH formation and evolution [TS] [IHl dZl dHj predict the 
existence of a large number of MBHBs emitting in the frequency range between 
~ 10^^ Hz and 10^^ Hz. PTAs can gain direct access to this population, addressing 
a number of unanswered questions in astrophysics (such as the assembly of galaxies 
and the physics of the relevant dynamical processes in galactic nuclei), by detecting 
gravitational radiation of two forms: (i) the stochastic GW background produced by 
the incoherent superposition of radiation from the whole cosmic population of MBHBs 
and (ii) GWs from individual sources that are sufficiently bright (and therefore massive 
and/or close) so that the gravitational signal stands above the root- mean-square 
(rms) value of the background. Both classes of signals are of great astrophysical 
and cosmological interest. In particular, the extraction of source parameters from the 
PTA observation may provide useful information about the system, helping in the 
identification of a putative electromagnetic counterpart. 

In this paper, we build on the results presented in [111 [20], review the features 
of the signal emitted by the population of MBHB and investigate how accurately the 
parameters of individual binary systems can be measured. The plan of the paper is as 
follow. In Section 2 we describe the GW signal relevant to PTAs, treating separately 
the unresolved background and the resolvable sources. In Section 3 we consider several 
representative MBHB population models for generating the GW signals. Results, in 
terms of signal characterization and detectability, are discussed in Section 4, and 
in Section 5 we report preliminary results regarding the parameter estimation of 
individual sources. We briefly summarized our main findings in Section 6. 

2. Gravitational wave signals for pulsar timing arrays 

A detailed derivation of the GW signal generated by a cosmological population of 
MBHBs is given in [TH] . In the following we will assume each source to be a circular 
binary defined by a chirp mass M = M^^^M^^^/iMi + Ma)^/^ {M2 < Mi are the 
masses of the two black holes) and a restframe emitted frequency fr = {l + z)f (/ is the 
observed frequency and z is the binary redshift), which is twice the orbital frequency. 
Consider now a MBHB population described by the distribution d^ N/dzdA4d\nfr, i.e. 
the comoving number of binaries emitting in a given logarithmic frequency interval 
with chirp mass and redshift in the range [M,M + dAi] and [z,z + dz], respectively. 
The characteristic amplitude he of the GW signal is given by: 




(1) 



where 



Hfr) 



8^2/3 ^5/3 
101/2 dL{z) 



ft 



■2/3 



(2) 



Gravitational Waves and Pulsar Timing 



3 



is the leading Newtonian quadrupole order contribution to the total GW emission, 
whose sky and polarisation averaged strain amplitude is given by ^T\, and d[^{z) is 
the luminosity distance to the source. 

2.1. The unresolved background 

In observations with PTAs, radio-pulsars are monitored weekly for periods of years. 
The relevant frequency band is therefore between 1/T - where T is the total 
observation time, assumed, unless otherwise specified, to be 5 years throughout this 
paper - and the Nyquist frequency 1 / {2 At) - where At is the time between two 
adjacent observations -, corresponding to 5 x 10~^ Hz - fewxlO"^ Hz. The frequency 
resolution bin is 1/T. Given a resolution frequency bin, we can introduce the concept 
of stochastic background. Whether the superposition of many deterministic signals 
should be effectively considered a stochastic signal depends on the frequency range 
and the observation time. Here we will consider the usual (simplified) criterion that 
the stochastic level of the signal is the amplitude at which the total number of sources 
contributing to the frequency bin of width 1/T centered at / is ^ 1 . 

The source distribution d^ N / dzdAidlnfr in equation ([T|) can be evaluated 
numerically from models of MBHB formation and evolution (see Section 3), and is 
typically a steep function of A4 (see figure 5 in jT9]). with very few massive binaries. 
As a consequence, at each frequency, the signal consists of the superposition of many 
contributions from relatively light binaries, overwhelmed by rare massive systems. A 
pragmatic way to evaluate the stochastic contribution to the signal is the following. We 
generate 1000 Monte-Carlo realization of the signal, computing numerically the mean 
differential distribution of sources with respect to redshift, chirp mass and frequency, 
i.e., d^ N / dzdAidf . If we consider an observed frequency /, and integrate the number 
of sources emitting in the frequency interval [/, /+1/T] in a given mass range [Ai, +oo[ 
and over all redshifts, for any given /, we can always find a M such that 



The integral above identifies the sources in the population that do not contribute 
to the stochastic background, i.e., those massive rare MBHBs which are likely to be 
individually resolved, because there is on average less then one of them per frequency 
bin. The result depends on the observational time, since the longer the observation, the 
narrower the frequency bin and, accordingly, the lower the level of signal contribution 
that can be considered stochastic. Broadly speaking the stochastic level is obtained 
by introducing in equation ([T|) an upper cut-off in the jM-integral at each frequency, 
according to the condition given by equation The spectrum of the observable 
induced residuals in the array is then given by 



2.2. Residuals from individually resolvable sources 

The angle-averaged optimal signal-to-noise ratio (SNR) at which a signal from a 
MBHB radiating at (GW) frequency ~ f can be detected using a single pulsar is 




(3) 



^ibkg(/) = K{f)/{2nf). 



(4) 



(5) 
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In the previous expression St^^sif) is the root-mean-square value of the noise level 
Stn at frequency /, including the instrument noise and the background contribution 
given by equation 21 (.) represents the average over the source position in the sky 
and orientation of the orbital plane, and (5igw(/) is the characteristic amplitude of the 
timing residual over the observation time T defined as: 



where the numerical pre-factor comes from the angle average of the amplitude of the 
signal and the y/jT term takes into account the build-up of coherent signal-to-noise 
ratio with the number of cycles. Equation ([5]) is appropriate to describe observations 
using a single pulsar; adding coherently the residuals from several pulsars yields an 
increase in SNR proportional to the square-root of the number of pulsars used in the 
observations. We will use the characteristic amplitude of the residuals (5ig„ to quantify 
the strength of a GW signal in PTA observations; Stgvj can be used to compute in a 
straightforward way the SNR, as a function of the noise level and number of pulsars in 
the array (all of which are quantities that do not depend on the astrophysical model), 
and therefore assess the probability of detection of sources in the context of a given 
MBHB assembly scenario. The statistics of individually resolvable sources is studied 
by running 1000 Monte-Carlo realisations of the whole population of MBHBs and by 
selecting only those sources whose characteristic timing residuals, given by equation 
([6|) , exceed the stochastic background level as defined by equations ^ and (jS]) in the 
previous section. 

3. Models of massive black hole binary populations 

The MBHB distribution (PN/dzdAidlnf,. is the only ingredient that is needed to 
compute both the background level and the statistics of the individual sources. We 
generate distributions of coalescing MBHBs using merging galaxy catalogs derived 
from the on-line Millennium Run database. The Millennium Simulation |22) covers 
a volume of (500//iioo)^ Mpc^ (here -ffioo is the value of the Hubble parameter 
today normalised to 100 km/s/Mpc) and is the ideal tool to construct a statistically 
representative distribution of massive low/medium-redshift sources (which were 
demonstrated to dominate the signal in the PTA window [IH])- As a first step we 
compile catalogues of galaxy mergers from the semi-analytical model of Bertone et al. 
[23] applied to the Millennium run. We need to associate to each merging galaxy in 
our catalogue a central MBH, according to some sensible prescription. The Bertone 
et al. catalogue contains many properties of the merging galaxies, including the bulge 
mass Mbuige, and the bulge rest frame magnitude My both of the progenitors and of 
the merger remnant. This is all we need in order to assign a MBH to each galaxy. 
The process is twofold. 

(i) We populate the coalescing galaxies with MBHs according to four different MBH- 
host prescriptions: 

• Mbh — Mbuigc in the version given by Tundo et al. ( 24|, "Tu" models); 

• Mbh — -Wbuige, with a redshift dependence in the version given by Mclure et 
al. ( [25], "Mc" models); 

• -^^BH — My as given by Lauer et al. ( [26!, "La" models); 

• -^^BH — cr as given by Tremaine et al. ( [27 , "Tr" models); 




(6) 
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To each merging system we assign MBH masses according to one of the selected 
models so that we have the masses of the two MBH progenitors, Mi and M2. For 
each prescription we also calculate the mass of the MBH remnant, M^, using the 
same relations. In all cases, the remnant mass is Mr > Mi + M2, consistent with 
the fact that MBHs are expected to grow predominantly by accretion during the 
merger process. 

(ii) For each MBH-host relation we consider three different accretion scenarios: 

• The masses of the coalescing MBHs are Mi and M2; i.e., accretion is triggered 
after the MBHB coalescence. We label this accretion mode as "NA" (no 
accretion) . 

• Accretion is triggered before the MBHB coalescence and only the more 
massive MBH (Mi) accretes mass; in this case the masses of the coalescing 
MBHs are aMi and M2, where a = {Mr - M2)/Mi - 1. We label this 
accretion mode as " SA" (single BH accretion) . 

• Accretion is triggered before the MBHB coalescence and both MBHs are 
allowed to accrete the same fractional amount of mass; in this case the masses 
of the coalescing MBHs are /3Mi and l3M2, where /3 = Mr /{Mi + M2) - 1. 
We label this accretion mode as "DA" (double BH accretion). 

Assigning a MBH to each galaxy, we obtain a list of coalescences (labelled by 
MBH masses and redshift). Each event in the list is then properly weighted over 
the observable volume shell at each redshift to obtain the distribution d^N/dA4dzdtr, 
which is finally converted into d^ N / dzdAddlnfr using the standard quadrupolc formula 
relation for the frequency shift dfr/dtr (equation (8) in [T9]). 

4. Results 

Single Monte-Carlo realizations of the typical signal, in terms of hc{f), are plotted in 
figure [1] for the four selected MBHB population models, assuming the "SA" accretion 
mode. The signal spectrum can be decomposed in a "confusion component" (the 
stochastic contribution), plus a handful of spikes due to particularly bright sources 
that are individually resolvable above the background level. In the following, we 
characterize both contributions. 

4-1. Background: characterization and observability 

The background level, computed following the prescription given in Section 2.1, 
deviates significantly from a simple f~^^^ power-law (see, e.g., 28 ) for / > 10^® 
Hz, and it is well fitted by a function of the form 



with ho = (1.46 ± 0.67) x IQ-^^, /q = 3.72j:^j2 x IQ-^Hz, and 7 = -1.08j:°:g!^. In 



general, the slope of the stochastic contribution to he starts to deviate from —2/3 at 
around 10~^ Hz, and becomes as steep as « —1.5. Most of the signal predicted by 
semi-analytical approaches (e.g. |3]) is actually not present because of the discrete 
nature of sources (see |19| for a detailed discussion). The errorbars are computed 
taking into account for uncertainties related to: (i) the galaxy merger rate at low 
redshift which determines the abundance of MBHBs, (ii) the MBH-host galaxy relation 




(7) 
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Figure 1. The characteristic amplitude of the GW signal from the population 
of MBHB systems. In each panel the thin lines identify the estimated level of 
the stochastic contribution assuming "SA" (solid line), "DA" (dashed line) and 
"NA" (dotted line) accretion modes. The total GW amplitude from a single 
Monte-Carlo realisation of the signal corresponding to the "SA" accretion mode 
is also shown as a thick solid line. 

assumed in populating galaxies with MBHs and (iii) the adopted accretion mode. The 
uncertainty range can be then compared to the sensitivities of ongoing and planned 
PTAs, as shown in Figure [2j The current sensitivity (as presented in [29]) is a factor 
of ^ 4 above the most optimistic estimate of the background. However, the complete 
Parkes PTA (PPTA) should provide a sensitivity improvement of about an order of 
magnitude in amplitude; several scenarios and/or parameter values within a given 
model could yield a detection at this level. The planned SKA will provide a major 
leap in sensitivity: a monitoring of 100 millisecond pulsars for 10 yrs at a precision 
level of i5trins ~ 50 ns will allow us to study the background spectral features in the 
frequency range 3 x 10^^ ^ / ^5 few x 10~®, enabling the complete characterization 
of the GW signal. 

4-2. Properties and statistics of resolvable binaries 

To quantify the statistics of the individual sources, we cast the results in terms of the 
cumulative number of resolvable sources as a function of the timing residuals: 

f°° dN 

^(^^8") = / ;77A7rT'^('54w) , (8) 

where the integral is restricted to the sources that produce residuals above the rms level 
of the stochastic background. Each Monte-Carlo realisation clearly yields a different 
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Figure 2. The sensitivity of PTAs to a GW stochastic background. The 
long-dashed thin lines represent the sensitivities of current and future timing 
experiments, as labeled in the figure. The thick lines depict the stochastic signal 
predicted by two representative MBH evolution models, see equation {Tjl. The 
shaded area marks the possible range of GW background level, according to the 
model uncertainties. 



value for N(5tgvj) ; the values shown in this section refer to the sample mean computed 
over the set of Monte-Carlo realisations. 

In figure |3] we show the mean total number of individual sources that exceed a 
given level of timing residual (as defined by equation (|5])), as a function of the timing 
residual. The left panel of figure [3] shows results where the accretion prescription is the 
same ("SA"), but the underlying MBH-host relation changes; assuming a sensitivity 
threshold of 30 ns, one expects to observe of the order of one source in the La-SA 
model, while there is only a probability « 5% for the Tr-SA model. The right panel 
shows that different accretion scenarios may cause a fiuctuation of a factor of ~ 3 
in the mean number of expected sources resolvable at a 10 ns level (e.g., between 
0.5 (Mc-NA) and 1.6 (Mc-DA)). In turn, the timing precision required for positive 
detection is in the range 5 — 50ns. 

5. Parameter estimation of resolvable systems 

5.1. Fisher information matrix formalism 

If a source is individually resolved, then the observations provide a means to measure 
its astrophysical parameters. Most of the signals detected by PTAs may be from 
circular binaries exhibiting negligible frequency drift during the observation time. For 
this model, the timing residuals ra{t; A) from each pulsar (labelled by a) depend on 
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Figure 3. Left panel: The effect of the MBH-host galaxy relation, assuming 
that accretion always takes place for a single black hole before merger ("SA" 
models), on the number of observable systems. The plot shows the number of 
total (thin lines) and resolvable (thick lines) sources 7V((5tgw) as a function of 
5tgw, see equation [S] Four different MBH merger scenarios are considered: Tu- 
SA (solid line), Tr-SA (long-dashed lines), Mc-SA (short-dashed lines) and La-SA 
(dotted lines), see Section 3 for a description of the models. Right panel: The effect 
of the MBH accretion model on the number of individually observable systems. 
The plot shows the number of resolvable sources only, N{Stg„) as a function of 



Sts 



As reference for the MBH-host galaxy relation, models "La" (thick lines) 



and "Tr" (thin lines) are considered. The line style is as follow: model La-SA 
and Tr-SA (solid lines), La-DA and Tr-DA (short-dashed lines), and La-NA and 
Tr-NA (long-dashed lines) . The duration of the observation is set to T = 5 yr . 



a 7-diinension parameter vector A {R,9, cfijip, t, f,^o}. Here, R ^ h/(2Trf) is the 
amplitude, 9 and (f> describe the source position in the sky, tp is the source polarization 
angle, t its inclination, / its (constant) frequency and $o its initial phase. Given an 
array of pulsar we can then compute the accuracy in the parameter estimation by 
computing the expected variance cr| = (r~^)^.^. of the parameter measurements, where 
r^fc is the Fisher information matrix for the observation in the array of a = 1, ...,N 
pulsars, and is given by: 



dra(t; A) 




dra{t] A) 



(9) 



Here, {x\y) denotes the noise weighted inner product of the two functions x and y - in 
our analysis we assume that all the pulsars have the same noise level, which is modeled 
as Gaussian, stationary noise - see e.g. [30]. With this definition, TinRjnR ~ P^, where 
p is the coherent PTA optimal signal-to-noise ratio. We normalized our calculation to 
a total SNR= 10 and assume a ten years observation for this exercise. The accuracy 
in the parameter estimation scales linearly with the SNR. 

5.2. Results for an array of pulsar 

We consider two different observation scenarios: (i) we consider a PTA consisting 
of N pulsars randomly and isotropically distributed in the sky varying N = 
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Figure 4. Left panel: Median errors in the source parameter estimation. Each 
triangle is obtained averaging over > 2.5 X 10* Monte Carlo generated sources. 
In each panel, solid lines (triangles) represent the median errors, and the thin 
dashed lines label the 25**^ and the 75"^ percentile in the error distributions. Left 
panel: errors as a function of the number of pulsars M, assuming a randomly 
distribution in the sky and a total SNR= 10 in the array. Right panel: Same as 
the left panel but now, in the a;-axis we vary the sky coverage of the array in the 
range [0,47r]rad, keeping A'^ = 100 and SNR= 10, linestyle as in the left panel. 



3, 5, 10, 20, 50, 100, 200, 500, 1000; (ii) we fix iV = 100, and we change the sky coverage 
of the PTA (defined as the minimum solid angle containing all the pulsars) Afijv = 
0.21,0.84, 1.84, 7r,27r,47r. For each pulsar distribution we generate > 25000 GW 
sources uniformly distributed in the sky with random initial phase and polarization, 
and inclination sampled according to a probability distribution = cos(t) in 

the interval [0, tt]. As a figure of merit of the performance of PTA we consider 
the median values of the parameter estimation accuracy as a function of N and 
Ar^AT. Results are shown in figure IH where we plot median and 25th and 75th 
percentile values of ajj . At least three pulsars are necessary to constrain the source 
parameters. Assuming an isotropic distribution of pulsars (left panel), at a fixed 
SNR, the parameter determination accuracy improves a lot with increasing N up 
to ~ 20. Adding further pulsars only improves the sky location of the source 
An = 27ry(shIM^A^)2^^^ihI^c^ (c'''^ is the 6* - correlation coefHcient [30]); its 
median value is well matched by the function Afi « Att/{p'^^/N). Assuming a plausible 
SKA configuration of 100 pulsar isotropically distributed in the sky, a source with a 
signal-to- noise ration of 10 can be located within an error box « 40 deg^; l and ip are 
determined within a 0.3rad confidence. The frequency is pinned down to sub-frequency 
bin resolution, in this case « 0.1 nHz. The right panel highlights the beneficial effect 
of having pulsar isotropically distributed in the sky, and not concentrated in a limited 
area. Although the errors on most of the parameters are insensitive to the array sky 
coverage as long as Arjjv > 1 rad, the sky position accuracy increases linearly with 
the sky coverage over the whole range [0, 47r]rad. Determining the source sky location 
is crucial to the task of finding a possible electromagnetic counterpart, which would 
have enormous scientific payoff. A uniform pulsar distribution in the array minimizes 
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also the probability of detecting a binary in a region of the sky where the PTA sky 
resolution is poor. 

6. Summary 

We have explored the GW signal generated by a cosmological population of circular 
MBHBs, and its detectability with ongoing and planned PTAs. The signal spectrum 
can be decomposed in a stochastic contribution, on top of which rare, bright and 
individually resolvable sources sit. The stochastic background is well described by a 
double power law function. The slope of the he spectrum starts to deviate from —2/3 
at around 10^® Hz, and becomes as steep as « —1.5. The current PTA sensitivity 
is a factor of ~ 4 above the most optimistic estimate of the background level. The 
complete Parkes PTA should guarantee a positive detection of the background for 
several selected scenario, while the planned SKA will allow to study its spectrum in 
the frequency range 3 x 10~^ ^ / ^ few x 10~^. A timing precision of 5 — 50ns 
(within the capabilities of the Parkes PTA and of the SKA) is required for detecting 
individual sources rising above the background. All the explored models predict at 
least 5-to-lO resolvable sources at a timing precision of Ins. We have also evaluated the 
accuracy with which the parameters of a source can be measured, using an analysis 
based on the computation of the Fisher information matrix. For a plausible SKA 
configuration, with 100 pulsars uniformly distributed across the celestial sphere, the 
source sky location can be determined within an errorbox of ~ 40deg^, assuming a 
total signal-to-noise ratio of 10 in the array. The detection of these systems, along with 
a putative electromagnetic counterpart, would provide invaluable information for the 
study of MBHB pairing and evolution at low redshift and their role in the formation 
and evolution of massive galaxies and clusters. 
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